A direct solution to the phonon Boltzmann equation 
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Abstract 

The frequency dependent phonon Boltzmann equation is transformed to an integral equation 
over the irreducible part of the Brillouin zone. Simultaneous diagonalization of the collision kernel 
of that equation and a symmetry crystal class operator allow to obtain a spectral representation 
of the lattice thermal conductivity valid at finite frequency. The method is applied to C, Si and 
Mg2Si to obtain the static and dynamical thermal conductivities. 
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The study of lattice heat transport in a crystal compound requires the knowledge of the 
phonon excitations as well as a model of transport which, in bulk systems, is conveniently 
taken to be the Boltzmann equation. Harmonic phonon spectrum are routinely obtained 
from ab-initio calculations based on the density functional theory, even for complicated 
compounds[l|. However to access transport properties like the thermal conductivity such 
calculations are not sufficient since it is necessary to describe the scattering of harmonic 
phonons by others phonons, impurities and crystal boundaries. Among those three scatter- 
ing processes the scattering by the others phonons is usually the more demanding since it 
originate in the anharmonic part of the total energy and therefore at least the third varia- 
tion of the energy with respect to atomic displacements is needed. Such calculations have 
however be shown to be feasible, either from density perturbation theory |2|, or using finite 
displacements [3|] . 

Once these scattering processes are calculated, and the related collision matrices constructed, 
the transport Boltzmann equation still remains to be solved. Among the few published re- 
sults for the ab-initio calculations of the thermal conductivity, all solve this equation itera- 
tively and consider the stationary case 4(] . It should be noticed that even in the simple case of 
silicon this scheme requires already several tens of iterations to get the desired accuracy. In 
addition there is no proof for the convergence of the iteration procedure. However in appli- 
cations there are many cases of interest where materials with a low thermal conductivity are 
needed and therefore where such an approach may be inadequate. It is for example the case 
in thermoelectricity where the lattice thermal conductivity should be as small as possible 
to increase the figure of merit ZT. In this context the study of time/frequency dependent 
lattice thermal conductivity is also important and is the subject of intense research^]. For 
a material to be a good thermoelectric it should have good electrical properties such that 
conductivity and thermopower but should also be a poor thermal conductor. These two re- 
quirements have been proved to be difficult to be achieved together in bulk materials. One 
way to circumvent the difficulty is to nano-striicturate the materialsy]. An other maybe 
to use finite frequency properties . It is known [7] that at a certain frequency the thermal 
conductivity start dropping rapidly what will increase the thermoelectric efficiency. The 
dropping frequency will be calculated here for the first time using ab-initio calculations. 
In this letter I present a direct non iterative solution to the Boltzmann equation applicable 
to the stationary and non stationary regime. This allow to calculate the static and dynam- 



2 



ical thermal conductivities and the accuracy of the method allows to study materials which 
are poor thermal conductors. A single parameter control the accuracy, mainly the number 
of points used to sample the first Brillouin zone. The paper is organized as follow. The 
Boltzmann equation is first reduced to an integral equation over the irreducible part of the 
Brillouin zone. Then the thermal conductivity is expressed in term of the collision operator 
defined during the reduction process. This operator is symmetric and, due to the reduction 
over the irreducible part of the Brillouin zone, is small enough to be diagonalized numeri- 
cally. Consequently a spectral representation is obtained for the thermal conductivity, valid 
at zero and finite frequency. Finally the method is applied to compounds having from very 
large to very low thermal conductivity. In each case the agreement with experiment is ex- 
cellent. 

The phonon Boltzmann equation is an integral equation over the first Brillouin zone which 



in its linearized version takes the form 
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The scattering of phonon appears in the right hand side of the equation through the term 
C{qp) for the collision processes, and with D{qp) for the decay processes. In the above 
equations ngp is the occupation function for a phonon of wave vector q in branch p. Vqp is 
the velocity and T = T(r, t) the temperature, nfp is the occupation function at equilibrium, 
and n^qp is the first order deviation from equilibrium, Uqp ~ n^p + n^p . It is possible to 
rearrange the scattering integral of ref 8| in order to make its relation to the lifetime of 
phonons calculated in js] more explicit jsj. 
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-^qq'Qb i^-.^^^ strength of the interaction in between the three phonons involved in the scat- 



tering [3[ and A a function which is zero unless its argument is a reciprocal lattice vector. 



in which case it take the value 1. 

We can then make the following ansatz for n^p, sinh (^-^^^^n^qj = fgp = 
Yla=i I dt' ^^-^ fqp{t — t'), where a is used to label the cartesian components of the vec- 
tor fqpit — t'). If equation [1] is transformed to Fourier space we obtain 
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The velocity is odd under inversion, v^qp = —Vqp, and it is easy to check that the collision 
matrix is even, fl'_gp _q,p/ = ^'qpq/pi, which means that f_qp{uj) = —fqp^u). Because the 
Brillouin zone contains q as well as —q it shows that the collision matrix is not unique 
and that it is indeed possible to make other choices flqp^q'p' such that '^q,p, ^'qp,q'p'fqp{^) = 
JZq'p' ^qp,q'p'fqp{^) ■ We choosc to work with Qqp^q'p' glvcu by 
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X [5{ujqipi — UJqp + COq^pJ + S{Uq'p' — Ugp — Wg^^pJ + 6{Uq>p> + Uqp — UJq^pJ] 

which can be obtained by a dummy change of variable q' — > —q' in the summation of equation 
|3l This matrix is clearly symmetric and can be shown to be positive definite using the same 
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method than in 
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In the following we denote by g a general point in the Brillouin zone, and by A; a point in 
the irreducible part of the Brillouin zone. R are rotations of the isogonal point group g of 
the crystal and \g\ denote the cardinal of that group. We denote by Qk the multiplicity for 
the branches of the star of k. 

In equation [3] if we restrict the velocity field to the irreducible part of the Brillouin zone the 
Boltzmann equation becomes 

TT- — - - \}hp,R'k'p' - tu}dk,R'k'Opp')-r-rjR'k'p'[^)- 

AksT^ sinh ^p, \9\ 

Under the rotations R the velocity transform like f^^^ = ^2/3 Rapv^p and it can be checked 
that the collision matrix is invariant, ^Rkp,Rk"p' = ^kp,k'p'- Therefore the Boltzmann equation 
written at point Rk show that f]^ik'p'{^) and J2i3 ^a^fRR'k'p'i^) fulfill the same equation. 
This gives f^^pi^) — ^2/3 -^apftpi^) + "^"(i^)) where u{u)) is any vector in the null space 
of flkp,R'k'p' — i'^^k^R'k'^pp' ■ We will see later that the null space of this operator does not 



contribute to the lattice thermal conductivity and therefore that we can choose u{uj) — 0. 

— * 

In other words fkp{^) transform like the velocity and can therefore be understood as being 
proportional to the phonon means free path. The Boltzmann equation can finally be reduced 
to an integral equation over the irreducible part of the Brillouin zone only, 
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In matrix notation this is written as \X) = —{fl — iu}P)\f{uj)) , with obvious definitions for 
VL and P, and 
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The operator P^ is working like the identity on vectors which transform like the velocity, 
Pk^kp — Vkp- is a collision matrix. It tells that when working in the irreducible part of 
the Brillouin zone, and considering a transition from vector k to /c', one should, obviously, 
also consider all the transitions to the different branches of the star of k' . Using the group 
properties of the set of rotation matrix, one can show that the matrices P^ and Q are 
symmetric. 

The energy flux through the lattice is given by JEif) = '^Yliqpf^qp^qp''^qp{^) therefore its 
Fourier transform is Je(u;) = —h{uj)^{ijj) with the thermal conductivity tensor given by 

qp sinh 2fcj^T ) iPi'p' 

In the second step we have used that the factor of fgp{uj) in the summand is just the drift 
term in the Boltzmann equation. As for the Boltzmann equation, the double integral over 
the Brillouin zone can be reduced to the irreducible part. We obtain 

with Xyfep,yfe'p'(Q;, ^) = Skk'Spp'^j^Ra-yRpy- This operator is diagonal in the kp space. Its 
value for the cartesian variables a and P depend on the symmetry class of the system and 



can easily be calculated using the great orthogonality theorem of groups theory. From its 
definition it is also clear that X(q;,/3) = X*(/3,q;) and that its purpose is to project out the 
components of the velocity and mean free path not involved it the a/3 component of the 
conductivity tensor. The operator Q{a, (3; u) = I{a, f3){Q — iuP) appearing in the thermal 
conductivity transform the same way as X, Q{a,P]u) = a; w), therefore using the 

symmetry of Cl and P we obtain the commutation relation [T{a, (3),Cl — iuP] = 0. This last 
identity implies the Onsager reciprocity relations at finite frequency, k^'^Iu) = /t^"(ci;) and 
therefore allow to obtain a more symmetric equation for the thermal conductivity, 

«:"''M = ^^(/M|(X(a,/3)+X(/3,a))(f2-za;P)|/H), 

because X(a, /3) + X(/3, a) is a symmetric matrix. It shows also that the null space of 
(X(a,/3) + X{l3,a)){Q — iuP) does not contribute to the lattice thermal conductivity. A 
vector which belongs to the null space of Q — iuP also belongs the null space of (X(a, (3) + 
X(/3, a)){Q — iuP). Therefore we can choose a solution to the Boltzmann equation which is 
orthogonal to ker(r2 — iuP) and transform like the velocities, = —{Cl — za;P)~^|X). 

Here ~ 1 is used to denote the Moore-Penrose inverse. The thermal conductivity can now 
be expressed as an average value over the known vector |X), 

K'^^ioo) = -^{X\{n-tujPr\l{a,^)+I{^,a))\X}. 

The matrices X(a, /3) +X(/3, a) and Q are symmetric and commutes. It is therefore possible 
to find a set of eigenvectors \er) such that Q\er) = ujr\er) and (X(a,/3) + X{(3,a))\er) = 
ir{<y, P)\er) ■ This gives finally a spectral representation for the dynamical thermal conduc- 
tivity, 

V ^ Ur — ico J u' — iuj 

r 

where papiuj') is a spectral density and the prime in the summation tells that the null space 
have to be excluded. 

The method has been applied to materials with high (diamond), medium (silicon) and low 
(magnesium silicide) thermal conductivity. Ah initio calculations are performed to obtain 
the interaction strength in between the phonons, F^^,^^, using the method in Q]. Therefore 
no adjustable parameters are used in the calculation. Al5xl5xl5x mesh is used to 
sample the Brillouin zone and the scattering by isotopes and surfaces have been included in 



the usual ways {4]. The results of the calculations are shown on figured^ for the static ther- 
mal conductivity. In each case the agreement with experiment is excellent. This has to be 
related to the use of the previous equation which only involves the diagonalization of small 
matrices and where symmetry has been used at best to reduce numerical uncertainties. 
The real (k^) and imaginary parts of the dynamical lattice thermal conductivity at 
room temperature are show on figure [T]d as a function of frequency. A rapid drop of 
is observed after some cut off frequency I/tq which also correspond to a maximum in Kj. 
Such rapid decrease for has already been obtained in silicon using molecular dynamics 



calculations 



and is known for long time in principle. When considering rapid time 



variations, the Fourier's law has to be modified to account for the finite time needed to 
establish a current, what may eventually lead to thermal waves at large frequency. More 



work is needed to obtain the phonon second sound of [12|, but fortunately the structure of 
the equations does not change and our procedure can be applied. Such work is in progress. 
The evolution of I/tq with temperature is shown in the inset of figure [1^. It is clear that 
at high enough temperature it becomes proportional to the temperature, as required for a 
relaxation time. It can also be seen on figure lb that the spread of Hi looks broader for 
silicon. This may indicate that a single relaxation time is not sufficient to account for the 
dynamics and would leads to a more complicated evolution of the temperature. This is 
indeed confirm by the spectral density Pa/B ploted in the inset of figure lb since there is a 
significant weight around 5 lO^Hz. 

To summarize, exploiting the symmetry of the system we have given a solution to the Boltz- 
mann equation which allow to compute the thermal conductivity from a spectral represen- 
tation. This way numerical errors are greatly reduced and therefore the study of material 
with low thermal conductivity is possible. In addition it is no more difficult to obtain the 
dynamical thermal conductivity which is calculated here for the first time. This allow for a 
quantitative estimate of the drop frequency of Kr- This can be of great interest in thermo- 
electric applications where k need to be reduced to increase ZT, but also for the industry 
of microprocessors. Considering their clock rate, heat transport at high frequency need to 
be understood. The drop frequency of the materials calculated here are very large. How- 
ever a decrease of Hr at much lower frequencies has been reported experimentally in alloys 
compounds and could lead to applications. 
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FIG. 1: Lattice thermal conductivities for diamond (red), Si (blue) and Mg2Si (green). The panel 
a) shows the calculated (full lines) and experimental (circles) static thermal conductivity as a 
function of temperature. Panel b) shows the real (full lines) and imaginary (dashed lines) parts 
of the dynamical lattice thermal conductivity. The inset in panel a) shows the inverse relaxation 
time To as a function of temperature. The inset in panel b) shows spectral densities. The full line 
represent while for the dashed line is just cx 6{uj — cOr)- The experimental data are taken 
from '^Ref[13] ''Ref|l4l '^RefQ ''RefQ- 
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